home *** CD-ROM | disk | FTP | other *** search
/ Over 1,000 Windows 95 Programs / Over 1000 Windows 95 Programs (Microforum) (Disc 1).iso / 1249 / secant.t < prev    next >
Encoding:
Text File  |  1997-04-18  |  2.3 KB  |  123 lines
  1. %
  2. % "secant.t" solve for multiple complex roots
  3. % using the secant method
  4. %
  5. %   Sample program for the T Interpreter by:
  6. %
  7. %   Stephen R. Schmitt
  8. %   962 Depot Road
  9. %   Boxborough, MA 01719
  10. %
  11.  
  12. const MAXIT : int := 20
  13. const ROOTS : int := 3
  14. type carray : array[ROOTS] of complex
  15.  
  16.  
  17. program
  18.  
  19.     var init, step : carray
  20.     var root, this, prev : complex
  21.     var i : int
  22.  
  23.     % you must set initial search points
  24.  
  25.     init[0].re := 1.0
  26.     init[0].im := 0.0
  27.     step[0].re :=-0.1
  28.     step[0].im := 0.0
  29.     
  30.     init[1].re := 0.0
  31.     init[1].im := 1.0
  32.     step[1].re := 0.1
  33.     step[1].im := 0.1
  34.  
  35.     init[2].re := 0.0
  36.     init[2].im :=-1.0
  37.     step[2].re :=-0.1
  38.     step[2].im :=-0.1
  39.  
  40.     for i := 0...ROOTS-1 do
  41.  
  42.         this.re := init[i].re
  43.         this.im := init[i].im
  44.         prev.re := init[i].re - step[i].re
  45.         prev.im := init[i].im - step[i].im
  46.  
  47.         search( root, this, prev )
  48.  
  49.         put "root =", cstr( root )
  50.  
  51.     end for
  52.         
  53. end program
  54.  
  55. %
  56. % evaluate function: fn(x) = x^3 - 2 x^2 + 4 x - 8
  57. %
  58. procedure fn( var dest, srce : complex )
  59.  
  60.     var temp3, temp2, temp1, temp : complex
  61.     
  62.     % x^3 term
  63.     cmov( temp3, srce )
  64.     cmul( temp3, srce )
  65.     cmul( temp3, srce )
  66.     
  67.     % x^2 term
  68.     cmov( temp2, srce )
  69.     cmul( temp2, srce )
  70.     temp2.re := temp2.re * 2
  71.     temp2.im := temp2.im * 2
  72.  
  73.     % x term
  74.     cmov( temp1, srce )
  75.     temp1.re := temp1.re * 4
  76.     temp1.im := temp1.im * 4
  77.  
  78.     cmov( temp, temp3 )
  79.     csub( temp, temp2 )
  80.     cadd( temp, temp1 )
  81.     temp.re := temp.re - 8.0
  82.  
  83.     cmov( dest, temp )
  84.  
  85. end procedure    
  86.  
  87. %
  88. % search for a root using the secant method
  89. %
  90. procedure search( var root, this, prev : complex )
  91.  
  92.     var dx, dfn, fn_this, fn_prev, next, delta : complex
  93.     var i : int
  94.  
  95.     for i := 0...MAXIT do
  96.  
  97.         cmov( dx, this )
  98.         csub( dx, prev )
  99.  
  100.         fn( fn_this, this )
  101.         fn( fn_prev, prev )
  102.         cmov( dfn, fn_this )
  103.         csub( dfn, fn_prev )
  104.  
  105.         if cabs( dfn ) < 1.0e-99 then
  106.             exit
  107.         end if
  108.  
  109.         cmov( delta, fn_this )
  110.         cmul( delta, dx )
  111.         cdiv( delta, dfn )
  112.  
  113.         cmov( next, this )
  114.         csub( next, delta )
  115.  
  116.         cmov( prev, this )
  117.         cmov( this, next )
  118.  
  119.     end for
  120.  
  121.     cmov( root, next )
  122.     
  123. end procedure